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We study the evolution of cosmological perturbations, using a hybrid approximation scheme which 
upgrades the weak-field limit of Einstein's field equations to account for post-Newtonian scalar and 
vector metric perturbations and for leading-order source terms of gravitational waves, while including 
also the first and second-order perturbative approximations. Our equations, which are derived in 
\f} • the Poisson gauge, provide a unified description of matter inhomogeneities in a Universe filled with 

' a pressureless and irrotational fluid and a cosmological constant, ranging from the linear to the 

highly non-linear regime. The derived expressions for scalar, vector and tensor modes may have 
£SJ ■ a wide range of cosmological applications, including secondary CMB anisotropy and polarization 

effects, cosmographic relations in a inhomogeneous Universe, gravitational lensing and the stochastic 
gravitational- wave backgrounds generated by non-linear cosmic structures. 
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> 1 I. INTRODUCTION 

7— I \ 

\^ , In recent years important results have been obtained in Cosmology from several observations of Cosmic Microwave 
£-■«. ■ Background (CMB) anisotropy and polarization, weak gravitational lensing effects and large-scale structure by means 
of galaxy redshift catalogs. The increasing precision that has been, and will be achieved by future experimental 
determinations requires comparable accuracy in the theoretical estimate of the several contributions to these effects. 
Accurate study of the evolution of cosmological metric perturbations is therefore crucial for understanding these 
contributions. Different kinds of techniques have been developed for this analysis, depending on the specific range 
of scales under consideration. For example, on scales well inside the Hubble horizon, but still much larger than the 
Schwarzschild radius of collapsing bodies, the study of gravitational instability of collisionless matter is performed 
■ using the Newtonian approximation. It consists of inserting in the line-element of a Friedmann-Robertson- Walker 
(FRW) background the lapse perturbation 2ipjy/c 2 , where ifN is related to the matter density fluctuation Sp via the 
C3 cosmological Poisson equation X7 2 tpN = ^Ga 2 5p. The dynamics of the system is then studied in Eulerian coordinates 
by accounting for the Newtonian mass and momentum conservation equations, owing to the fact that the peculiar 
matter flow v never becomes relativistic [1]. The Newtonian limit, according to which the gravitational field ip is 
always much less than the square of the speed of light, c 2 , can be improved by a post-Newtonian (PN) approach 
to account for the moderately strong gravitational fields generated during collapse. In this case, by considering the 
expansion of the general relativistic equations in inverse powers of the speed of light, it is possible to neglect terms 
of order (v/c) 4 in the equations of motions, i.e. to perform a first post-Newtonian (IPN) approximation, which in 
Eulerian coordinates accounts for non vanishing shift components and for an extra perturbation term — 1(pjq j(? in the 
spatial part of the line-element. Calculations using higher and higher orders of 1/c would generally lead to a more 
accurate description of the system, e.g. accounting for the generation of gravitational waves, and possibly allow for an 
extension of the range of scales to which the formalism can be applied. A PN approach to cosmological perturbations 
has been followed in Refs. [2-4], using Eulerian coordinates, while Ref. [5], uses Lagrangian coordinates. 

On the other hand, the first-order perturbations for non-relativistic matter, obtained with the Newtonian treatment 
coincide with the results of linear general relativistic perturbation theory in the so-called longitudinal gauge [6] . The 
relativistic linear perturbative approach is the one applied to the study of the small inhomogeneities giving rise to 
large-scale temperature anisotropics of the CMB. However, on small and intermediate scales, linear theory is no 
longer accurate and a general fully relativistic second-order perturbative technique is required. In fact, second-order 
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metric perturbations determine new contributions to the CMB temperature anisotropy [7, 8] and polarization [9]. 
In particular, second-order scalar, vector and tensor metric perturbations produce secondary anisotropies in the 
temperature and polarization of the CMB which are in competition with other non- linear effects, such as that due 
to weak gravitational lensing produced by matter inhomogeneities, which induces the transformation of E-mode into 
B-mode polarization [10]. 

Moreover, accounting for second-order effects helps to follow the gravitational instability on a longer time scale 
and to include new non-linear and non-local phenomena. The pioneering work in this field is due to Tomita [11] who 
performed a synchronous-gauge calculation of the second-order terms produced by the mildly non-linear evolution of 
scalar perturbations in the Einstein-de Sitter Universe. An equivalent result, was obtained with a different technique 
in Ref. [12]. The inclusion of vectors and tensor modes at the linear level, acting as further seeds for the origin 
of second-order perturbations of any kind was first considered in Ref. [13]; in Ref. [14] the evolution of relativistic 
perturbations in the Einstein-de Sitter cosmological model was considered and second-order effects were included in 
two different settings: the synchronous and the Poisson gauge. 

As we have stressed, the evolution of cosmological perturbations away from the linear level is rich of several effects 
as, in particular, mode-mixing which not only implies that different Fourier modes influence each other, but also that 
density perturbations act as a source for curl vector modes and gravitational waves. 

The aim of the present paper is therefore to obtain a unified treatment able to follow the evolution of cosmological 
perturbations from the linear to the highly non-linear regime. As we will show hereafter, this goal is indeed possible 
on scales much larger than the Schwarzschild radius of collapsing bodies, by means of a "hybrid approximation" 
of Einstein's field equations, which mixes post-Newtonian (PN) and second-order perturbative techniques to deal 
with the perturbations of matter and geometry. In our study we adopt the Poisson gauge [15], which, being the 
closest to the Eulerian Newtonian gauge, allows a simple physical interpretation of the various perturbation modes. 
We derive a new set of equations which holds on all cosmologically relevant scales and allows to describe matter 
inhomogeneities during all the different stages of their evolution. The new approach gives a more accurate description 
of the metric perturbations generated by non-linear structures than the second-order perturbation theory, which 
can only account for small deviations from the linear regime. For example, on small scales our set of equations can 
be used to provide a PN description of metric perturbations generated by highly non-linear structures within dark 
matter halos, while describing their motion by means of the standard Newtonian hydrodynamical equations. On 
large scales our equations converge to the first and second-order perturbative equations as obtained in Ref. [14] (see 
also Refs. [16, 17]), which implies that they can be applied to every kind of cosmological sources. 

The plan of the paper is as follows. In Section II we obtain metric perturbations in the Poisson gauge according to 
our hybrid approximation. In Section III we obtain the source terms for scalar, vector and tensor metric perturbations 
and the evolution equations for the matter density and peculiar velocity. Section IV is devoted to a comparison of 
our approach with known approximation schemes, such as the standard Newtonian one, the linear and the second- 
order relativistic perturbative approaches, in the appropriate regimes of applicability. Moreover, we derive a PN 
approximation to describe the highly non-linear regime on small scales. In Section V we show that the PN expressions 
for vector and tensor modes actually hold on all cosmological scales. In Section VI we sketch some cosmological 
applications of our approach and we draw our main conclusions. An Appendix is devoted to the solution of the 
inhomogeneous gravitational-wave equation. 

II. SCALAR, VECTOR AND TENSOR METRIC PERTURBATION MODES 

Adopting the conformal time 77, the perturbed line-element around a spatially flat FRW background in the Poisson 
gauge [6, 14] takes the form 

ds 2 = a 2 {ri) [-(1 + 24))drf - 2V a dr]dx a + ((1 - 2tp) 8 a p + h afj ) dx a dx (3 ] . (1) 

In Eq. (1) the metric includes perturbative terms of any order around the FRW background. In this gauge, V a 
are pure vectors, i.e. they have vanishing spatial divergence, d a V a = 0, while h a p are traceless and transverse, 
i.e. pure tensor modes, h a a = d a h a p — 0. Here and in what follows spatial indices are raised by the Kronecker 
symbol S a p. Unless otherwise stated, we use units with c = 1. From the results of the post-Newtonian theory 
[5, 18-20], we deduce that vector and tensor metric modes, to the leading order in powers of 1/c, are respectively of 
C(l/c 3 ) and 0(l/c 4 ), since to the lowest order it is well known that the line-element (1) assumes the weak-field form 
ds 2 = a 2 [— (1 + 2</>i) di] 2 + (1 — 2^1) S aj 3dx a dx t3 }, where the scalars (f>i and ^1 are both of C(l/c 2 ) and (f>i = ipi = <p 
[!]• 
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Let us now write Einstein's equations G l j = n 2 T l j in the perturbed form 

f'G'j + tfCPj = k 2 ( (0 »r 3 + tfT*,- ) , (2) 

where k 2 = 87rG/c 4 and ^G l j = n 2 ^T l j reduce to the background Fricdmann equations. Hereafter, we assume 
that the Universe is filled with a cosmological constant A and a pressureless fluid - made of Cold Dark Matter (CDM) 
plus baryons - whose stress energy-tensor reads T % j = pu l Uj (u l Uj — —1). In this case the Friedmann equations read 
ZH 2 = a 2 (8irGp + A) and p ' = -3Hp, where primes indicate differentiation with respect to 77, TL = a' /a and p is the 
mean matter energy density. 

Since the metric (1) can be expressed in an invariant form both in the PN and in the second-order perturbative 
approximations, we introduce a hybrid formalism that consists in evaluating Einstein's field equations up to the correct 
order in powers of 1/c, while including some "hybrid" correction terms, which cannot be completely absorbed in a 
"rigid" PN approximation, but are required for consistency with a second-order general relativistic approximation. 

Namely, writing ip = -01 +"02 and = 0i + 2 (where tp 2 and 2 in principle contain all powers of 1/c), and replacing 
the metric (1) in Eq. (2), the (0-0) and (0-a) components of the perturbed Einstein equations take the form 

6G% = ~ [-6H (H<j)i + 0i) + 2V 2 V + 3(9^1 ^V-i + 80iV 2 0i + Si] = k 2 5T° , (3) 
where Si = 12W 2 2 + 3^i 2 + 12W0iV>i - 127^1^1, and 



SG\ = ~ (nd a ^ + d a ^ + hl 2 V a + S 2q ) = k 2 5T\ , (4) 

where S 2a = -47Y0i<9 Q 0i - 20! d a ip[ - ip[d a 0i + 20i<9 Q 0i + 2ip 1 d a ip' 1 . 
The traceless part of the (a-/?) Einstein's equations 5G a p = k 2 ST a p reads 

| v 2 (0 - v>) - 1 0iV 2 0r - 1 0! v 2 0! + 1 ^v 2 0! - 1 erfadvfa + \ ^0i^Vi - 2^v>i^V>i s a p 

- 2d a dp(<p - V) + 8i/)id a dpil; 1 + 6d a Vi3/#i - 2<9 Q 0i<9 / #i + 2d a (j> 1 dp<i> 1 - 2d a tp 1 dpc/) 1 + 40i<9 Q <9 (3 1 

- A^dp^ + d a (2HVp + Vp) + dp(2HV a + V a ) + h" a (i + 2Hh' a p - W 2 h a p = 2a 2 n 2 (5T a p - ±6T\5 a p) , 

(5) 

while its trace becomes 

2V 2 (0 - V) + 6W0'i + 6(H 2 + 2W')0i + 6< + lZHifa - 2V fadvfa 

- 40iV 2 0i - 3a„^ia"Vi - sv-iv 2 ^ + 4i/-iV 2 0i - 2<%0i<9 i >i = aVar",, . (6) 

The components of the perturbed stress-energy tensor 5T l j will be calculated later to the correct order in powers of 
1/c. Note that Eqs. (3), (5) and (6) are evaluated up to C(l/c 4 ), while Eq. (4) is evaluated up to C(l/c 3 ). However, 
the terms Si and S 2a are at least of 0(l/c 6 ) and C(l/c 5 ), respectively, and come out from our hybrid scheme, 
which mixes PN and second-order perturbative approaches. Moreover, since our purpose is to calculate the source 
of gravitational waves to the leading order in powers of 1/c, in Eq. (5) we do not take into account contributions of 
order higher than 0(l/c 4 ), though retaining time derivatives of h a p. 

Taking the divergence of Eq. (4), to solve for the combination H(f>i + ip[, and replacing it in Eq. (3), we obtain 

V 2 (H0x + v4) = -^d v 8T\ - d v S 2v , (7) 

and 

V 2 V 2 ^ - V 2 V 2 (^! + fo) - [v 2 <5T% + m{d v 8T° v )\ - md v S 2v - ^V 2 Si - V 2 Qd^i^i + 4</>iV 2 V>i) . 

(8) 

The pure vector part V a can be isolated by replacing V 2 (Hcfii + ip[ ) in Eq. (4) , where we now neglect the term S 2a 
since it is at least of C(l/c 5 ) 

V 2 V 2 V Q = 2aV {d a d u 5T° v - \7 2 ST° a ) . (9) 
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Finally, applying the operator d^d a to Eq. (5), we can solve for the combination <f> — ijj and write 

v 2 v 2 (^ - v) = -\a 2 n 2 d (i d a (sT a - ±5T\5 a ^j + ^d a d f3 {d^idPfa - dp^d^ - d^d^ - d^d"^) 

+ V 2 Qft^ifl"^ - 29^1^01 + ^faFfa + 4^i V 2 V>i + 20!V 2 0! - 2^i V 2 ^ . (10) 

Replacing the latter expression in Eq. (5), together with the expression for the vector mode V a obtained by taking 
the divergence of Eq. (5), we find 

V 2 V 2 (ti la p + 2Hh' a - V 2 h a p ) = 2K 2 a 2 [V 2 {V 2 n a ti - d a d v H v - d f3 d v n a v ) 

+ \ (v 2 d»d u n\ s a p + d a d^d v ii\)] , (11) 

where we have defined the traceless tensor 

fl a p = ST a p - \5T%5 a D - (\i>iV 2 ^ + \d v c$> x d v i> x - ^iVVi - ^iV 2 ^ - \d v - d'^d^) 6 a (j 

o K ffl \ o o o o o J 

- -FT (^i^^Vi + 3^1^! - SP<hdp^h. + d a 4>idp<t>i - d^dpfa + 2cj> 1 d a d <j> 1 - 2^ 1 d a d p 4> 1 ) . (12) 
k a 

The form of Eq. (9) and Eq. (11) allows to directly check that vector sources are transverse while tensor sources are 
doubly transverse and traceless. 

Actually, there is a very simple way of solving the perturbed Einstein equations SG° a — K 2 a 2 6T° a and SG a p = 
n 2 a 2 ST a p with respect to vectors V a and tensors h a p, respectively. In fact, after retaining only the metric terms 
which appear linearly on the LHS of Eq. (4) and Eq. (5) [30] and defining as lZ a p the consequently obtained RHS, it 
is possible to apply to both sides the correct combinations of the direction-independent projection operator [21] 

r a f s = 8 a -(V 2 y 1 d a d (13) 

and automatically obtain the vectors 

\7 2 V a = -2a 2 n 2 V\ (ST° V - -^S 2l ) , (14) 

and tensors 

7 2> a . _ o.,2„2 ' — - 1 



h" a p + 2Hh' a p - V 2 h a p = 2/sV \V a v V»p - -V a pV^j K\ . (15) 

After applying twice the Laplacian operator to Eq. (15) and neglecting as before the term <S 2 i> in Eq. (14), we recover 
Eq. (9) and Eq. (11). 

III. THE STRESS-ENERGY TENSOR AND THE SOURCE OF METRIC PERTURBATIONS 

For the purpose of calculating the components of the perturbed stress-energy tensor 8T l j to the correct order in 
powers of 1/c, it is convenient to restore the speed of light c in the time coordinate dx° = cdrj. From Eq. (1) we 
obtain the four- velocity u % = dx % jds 



1 / „ , v 



l -2{ 2 ^-* 



+ 0[-t), (16) 



va -u° + o(\), (17) 



c \ c 5 



where v 2 — v v v v and v a = dx a /dr] is the coordinate three-velocity with respect to the FRW background. 
The total energy-momentum tensor for our fluid of dust and cosmological constant reads 

T\ = {0) T k + ST l k = [(p A + p)c 2 + Pb ] gu^u 3 + Pb 5\ , (18) 
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where ^T l k is the background stress-energy tensor and p = p + 5p is the total mass-density. The background density 
Pb = p + Pa includes the contribution from the cosmological constant, p\ = (Ac 2 )/(87rG), while the background 
pressure Pb — Pa — —pAC 2 is only due to the latter. 

Turning to the components of the perturbed energy- momentum tensor, in terms of the coefficients 0iand xjji of the 
metric (1) and up to the correct orders in powers of 1/c, we find 



ST% = T° - (°)T% = -cHp -v 2 p + (Jz) , 



(19) 



5T\ - T\ - (°)T° Q = v aP c (1 - 2fc - 2Vi) + — pv 2 + O ' ' 



(20) 



grpa = T a Q _ (OJya = _ y a _ ^ 2 + ( J_ 

c V c 6 



(21) 



«^ = - C)^ = v>p (l - 20! - 2^i + J) + O 



(22) 



In the hybrid equations which we are about to derive we will keep step by step only the 5T l j components we need to 
let our set of equations hold in the first, second perturbative order and PN regimes. 
By substituting 5T a /) in Eq. (10), we obtain up to 0(l/c 2 ) 



V 2 V 2 (0i -Vi) = o, 



(23) 



and we can safely assume 0i = V>i = <P- 

This allows to further simplify Eqs. (6)-(10) and obtain our final set of hybrid equations for cosmological perturba- 
tions, namely 



a 2 n 2 c 2 



V 2 (7^ + ^') = s— 0"(/w„), 



(24) 



v 2 v 2 v> = — 



n 



3 ../; 



(25) 



V 2 V 2 y« = 2a 2 k 2 [d a d v (cpv v ) - V 2 (cpv a )] 



(26) 



V 2 V 2 (0 -i/>) = -^aVd'U ^pw"^ - Jp« 2 <5%) + V 2 (5^3^ + 4pVV) - 3d„d» (d^d"<p) , (27) 



V 2 V 2 0=- 



a 2 K 2 



3wa" (p« w (i - 4^)) + 3#uo> - v 2 (c 2 5p + 2p« 2 ) + -v 2 (a^a» 



3d v d"(d ll tpff'tp) + 3^V (^ipdvip - Vd„^ - v 2 (e-^V + 2 ) , (28) 



V 2 V 2 (± h " a p + ™h' a fj - v 2 h a f ^j = 2kV [v 2 (v 2 n a p - d a d„n v - d p d v n a v ) 



+ ^ (V 2 cPd„ft% S a + d a d^d v Tl\)] , (29) 
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where the traceless tensor lZ a /3, Eq. (12), now reads 



while the trace part of the (a-f3) component of Eq. (2) becomes 

2V 2 (</> - V) + 18-^' + 6 (2^ + ^-]p + %tp" - 79^a> - 8pVV = K 2 a 2 p« 2 • (31) 

Using Eqs. (19)-(22), and the expression (j>\ = ipi = tp, we can write the stress-energy tensor conservation equations 
T-.j = in a form that will give us the equations for our pressureless fluid in the first, second-order and PN regimes, 
respectively. 

More specifically, the energy conservation equation reads 

p' + 3Hp + d v {pv v ) --ptp' + \[(pv 2 Y + d v (v v pv 2 ) + AHpv 2 ] - 2pv v d v tp --ptp' 2 - -pW = , (32) 
c c z c c 

while the momentum conservation equation reads 

(pv a )' + d v (pv"v a ) + ATLv a p + pc 2 d a tp - 16Hpv a ip + \d v (pv 2 v u v a ) - 2pc 2 tpd a tp - HcpV a + pcv u d a V v 

- Ad u (pv"v a ip) + ^Hpv 2 v a - Ap'v a <p - Apv' a <p - 6pv a tp' + (v a pv 2 )' - 2pv a v u d u (p + 2v 2 pd a ip + pc 2 d a cj) 2 = 

(33) 

Eqs. (25), (26) and (28)-(29), together with Eqs. (24), (32) and (33), are the main result of our paper and represent 
a new set of equations which allow to describe the evolution of metric perturbation from the linear to the strongly 
non-linear stage in terms of the gravitational field ip, the matter density p and the peculiar velocity v a . 

IV. LIMITING FORMS OF THE HYBRID APPROXIMATION IN DIFFERENT REGIMES 

We now show how our approach accounts for known approximation schemes in different regimes. 

A. The linear perturbative regime 

Linearly perturbing Eqs. (24), (25) and (28) with respect to the FRW background we deduce that (f> and ip coincide 
and we obtain the linear scalar potential <p in terms of first-order density and velocity fluctuations, 

\7 2 \7 2 p = 4a 2 irG[\7 2 6p~3Hpd u v„] , (34) 

V 2 {TLp + tp') = -Aa 2 TiGpd v v v . (35) 
Moreover, linearizing Eqs. (26), (29) and (31), using the linearized expressions for ST l j we obtain 

ip" + iTiip 1 + (2H' + H 2 )ip = , (36) 

V 2 V 2 V a = (37) 

b!' (i a + 2Hti p a - V 2 h a p = . (38) 

Perturbing to first order Eqs. (32) and (33), we recover also the linear continuity and momentum equations which 
read respectively: 

bp' + ZHdp + pd v v u -Zpp' = 0, (39) 

v' a +Hv a +d a tp = 0. (40) 

In other words, we obtain all the results of linear perturbation theory (see e.g. Ref. [22, 23]), if we interpret tp as the 
linear scalar potential. 
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B. The second-order perturbative regime 

On the other hand, selecting only the growing-mode solution of Eq. (36) and perturbing up to second order Eqs. (25)- 
(29), in the limit of a pressureless and irrotational fluid, we recover all the results of second-order perturbation theory 
[9, 14, 24]. 

More specifically, the first-order vector metric perturbations vanish, while the linear tensor metric perturbations 
are negligible for every kind of cosmological sources, thus we can safely neglect terms which can be expressed as 
products of first-order scalar and tensor metric perturbations. Writing ip(x,r]) = ipo(~x.)g(ri), where (po is the peculiar 
gravitational potential linearly extrapolated to the present time and g = D + /a is the so-called growth-suppression 
factor, where D + (jj) is the linear growing-mode of density fluctuations in the Newtonian limit, and using the results 
of the previous subsection we obtain 



V 2 VVs = -AnGa 2 [3H {pd v v sv + 8pd"v v + v v d v Sp - 4pd l, (<pv l/ )) - 3pd v (ip'v„) - W 2 {S Ps + pv 2 )] 
'3 



v 2 ( \d v ipd v v + 4y>vV ) , (4i) 



V 2 V 2 c/> s = -AirGa 2 



m {pd v v sv + 8pd v v v + v v d u 5p - ipd v (ipv u )) - 3pd v {ip'v v ) - V 2 {5 Ps + 2pv 2 ) 



+ 3pd»d v {v v v tl )\ + ^v 2 (d^d» - 3d^(d^d^) , (42) 

V 2 V 2 y Q = \&TiGa 2 d u (v„d a 5p - v a djp) , (43) 
V 2 V 2 (h" a p + 27ih' a fj - V 2 h a fj ) = 16 7 rGa 2 [V 2 (V 2 ^ a /3 - d a d v Tl v (i - d (i d v 1Z a v ) 

+ \ {v 2 d^d v n\ s a p + d a d^d u n\)] , (44) 

where R a p has the same analytic form of Eq. (30), that is 

K% = p (v a v p - l -v 2 5 a ^j - (d a <p d p <p - \d v <p d vV 5 a f^j - ^-i^ (v> d a d^ - ±<p V 2 ^ 5 a ^j , (45) 

The subscript s indicates quantities evaluated at the second perturbative order, vf is the velocity dx a /drj perturbed at 
the second-order and is related to the second-order spatial part v" 2 ^/a of the 4- velocity by the relation vfy = vf — <pv a 
[31]. 

We can also find the equations that describe the evolution of Sp s and vf by perturbing up to second order Eqs. (32)- 
(33) and taking the divergence of the latter. In this way we recover also the second-order energy continuity equation 

Sp' s + 3H6p s + pd v v v s + 5pd v v v + v v d v 5p + Hpv 2 + p{v 2 )' - 36p<p' - 3pip' s - Qptptp 1 - 2pv"d^ = , (46) 

and the divergence of the second-order momentum conservation equation 

Hpd a v sa + pd a v' sa + p V 2 cf) s + d a [4H6pv a - 4Hpifv a + 5p'v a + Spv' a - 4pLpv' a - §p(p'v a + Spd a ip 

- 2p<pd a <p + pd v {v v v a )] = . (47) 

C. The Newtonian approximation 

From Eqs. (23) and (25) up to £>(l/c 2 ) we deduce 

V 2 ^ = V 2 0! = W = ^S^Sp , (48) 

and writing ip = <pn/c 2 , we recover the Poisson equation V 2 (fiN — 4nGa 2 Sp, where the subscript TV stands for 
Newtonian. 
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Analogously, to leading order in 1/c, Eqs. (32) and (33) respectively become the usual continuity and Euler equations 
of Newtonian cosmology which apply in the limit of weak fields and non-relativistic velocities [6] 

f/ + 3Hp + d v (pv v )=0, (49) 

v' a + Ttv a + v v d v v a = -d a ip N . (50) 

The latter equation was obtained taking Eq. (33) up to O(l/c ) and inserting Eq. (49). 
In the linear limit Eqs. (48)- (50) become 

VV = 4nGa 2 6p , (51) 
Sp' + mbp + pd v v v = , (52) 

v' a + Hv a = -d a (p N . (53) 

As we can observe, the equations which characterize the linearized Newtonian theory, differ from the linearized 
relativistic ones. In particular, while the momentum conservation Eqs. (40) and (53) are identical, the linear energy 
density conservation Eq. (39) differs from the Newtonian one, Eq. (52), by the extra term —3pip' which does not vanish, 
even for the pure growing- mode solution of Eq. (36), owing to the presence of a cosmological constant contribution 
to the FRW background. 

Moreover, Eq. (34) represents the linear relativistic generalization of the Poisson equation, since it includes the 
contribution of the so-called longitudinal momentum density <pf (d v <ff = —4:a 2 irGpv„) which acts as a source term 
for the linear potential (p. Thus, the Poisson gauge gives the relativistic cosmological generalization of Newtonian 
gravity [6]. 

D. The highly non-linear regime in the PN approximation 

Finally, we consider the case of cosmic structures, in the highly non-linear regime, whose size is much larger than 
their Schwarzschild radius (in order to avoid non-Newtonian terms in the expressions of the sources) . 

Our sources can generate vector and tensor metric perturbations by mode-mixing in the non-linear regime. In 
particular, this mechanism applies to dark-matter halos around galaxies and galaxy clusters or, more specifically, to 
the highly condensed substructures by which these halos are characterized. 

We obtain the continuity and momentum equations up to 0(l/c 2 ), the equation describing the evolution of the 
(0-0) component of the metric (1) up to C(l/c 4 ), and the equation for the vector modes V a up to C(l/c 3 ), i.e. their 
1PN approximation. Moreover, we describe the scalar mode of the (a-(3) component of the metric (1) up to C(l/c 4 ), 
i.e. we consider its second post-Newtonian (2PN) approximation, while we obtain the leading-order terms in powers 
of 1/c for the source of the tensor modes h a p. 

Eqs. (25)- (29) in this limit become 

V 2 W = ^-V 2 Sp+ [VW) - 3Hd"(v v p)] lv 2 {^ip N d v ip N + 4^VViv) , (54) 

V 2 V 2 = ^-V 2 Sp + [2\7 2 (pv 2 ) - md v {v v p) - W x d v {pv v v ll )] + ^/ 2 {d v y N d v y N ) 



- ^d v & i (dr<p N ff'(pN) , (55) 

V 2 V 2 V4 = [d a r{v v p) V 2 (v aP )} , (56) 

V 2 V 2 (L h " a + ™h'% - V 2 h a ^j = [V 2 (V 2 TO - d^TVfs d p d»K%) 

+ \ {v 2 d»d„n\ s a p + d a dpd»d v K\)} , (57) 
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where the post-Newtonian limit of the traceless tensor lZ a 8 is 



, = ,>[v a vp-^5 a p 



1 



47rGa 2 



d a ip N dpip N - ^d"tp N d u <p N 5 a p 



1 



27rGa 2 



AttGu 2 



if N SpS a /3 

(58) 



Since, in order to compute the metric coefficients up to the PN approximation, we only need the terms in 8T % j which 
satisfy the Newtonian equations of motions, in Eq (58) we have inserted the Poisson equation. The 1PN extensions 
of the Newtonian continuity and Euler equations respectively read 



1 



p' + 3H P + d v {pv v ) + - [(pv 2 )' + d v {v v pv 2 ) + AHpv 2 - 2pv u d v ip N ] = , 



(59) 



p (v' a + Hv a + v v d v v a + d a (p N ) + — 



-4p'v a ip N - Apv' a tp N - 6pv a tp' N + (v a pv 2 )' - 2pv a v v d v Lp N + 2v 2 pd a (p N 



+ pd a 4>PN - 2p<PNd a (pN - HpV a + pv^daV^ - l6Tipv a tp N - Ady (pv u v a ip N ) + d v (pv 2 v u v a ) + ATLpv 2 v a 



0, 
(60) 



where 0pjv is given by the 1PN part of Eq. (55). It can be worth noting that the sources of the metric coefficients 
involve only quantities of Newtonian origin, i.e. they do not contain terms defined in higher-order approximations. 

To conclude this subsection, let us stress that all the PN expressions derived here are new, as they are derived in a 
different gauge than the usual post-Newtonian [2, 3, 18, 25], or synchronous and comoving one [5]. 



V. VECTOR AND TENSOR MODES 



It can be worth to observe that, in the linear limit, Eq. (31) becomes 

iffit + iW N + (2H' + H 2 )^ N = . 



(61) 



This result is extremely important since it implies that the Newtonian potential ifN and the linear potential ip evolve 
in the same way with time. Eq. (61) can be also obtained by mixing together the Newtonian continuity, Eulcr and 
Poisson equations perturbed at first order. This means that, starting from the same initial conditions, i.e. from the 
same primordial potential as given e.g. by inflation, the two linear potentials (Pn/c 2 and <p will assume the same values 
in each point and at each time. In other words, Eq. (61) implies that, in the case of first-order matter perturbations, it 
is sufficient to use Newtonian gravity on all scales, provided that we define a "Newtonian" linear density perturbation 
dpN via the Poisson equation applied to the linear relativistic potential ip, even if Sp^ differs from the relativistic 
density dp, as given by Eq. (34). The previous considerations allow to conclude that, for pure growing- mode solutions 
of Eq. (61), in the case of an irrotational and pressureless fluid, Eqs. (56)-(58) apply to all cosmologically relevant 
scales, i.e. from super-horizon to the smallest ones, even if the density p, the velocity v a and the potential ip are 
required to follow the usual Newtonian hydrodynamical equations. In the equations that follow, therefore, we will 
drop the subscript N on the various quantities and write 



VV = 4irGa 2 5p 



(62) 



P ' + mp + d v (pv v ) = o, 

v' a + TLv a + v v d v v a = -d a ip . 



Thus, for the vector modes we have 



and, for the tensor modes, 



WirGa 2 



(63) 
(64) 

(65) 



[V 2 (V 2 ft a /3 - d a d v Tl v fi - d d v H a v ) 



(66) 
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where 

K a p = p (v a vp - \v 2 6 a ^j - (&*<pdw - ±Br<pd v ip& a ^ - (vd a dp<p - ^-cpSpS^p 



(67) 



which represent the most important results of our paper, since these equations imply that, in the case of matter 
perturbations, the Newtonian description of the sources of vector and tensor metric fluctuations can take into account 
all the effects of the relativistic second-order perturbation theory. 

It is important to stress that the third term on the last line of Eq. (67) does not contribute to the source of 
gravitational waves since it vanishes after applying the projection operation in Eq. (13); thus we are allowed to drop 
it and define as effective source of the gravitational wave h a p the traceless tensor 



K a cS g = p [v a vp- -v 2 6% )+ —— d a V d^ - -8 "tp d vV 8 a p . (68) 



Actually from a post-Newtonian point of view, these equations hold true also for a pressureless fluid with a vorticity 
contribution to the peculiar velocity v a , but the reader should not be surprised if curl terms can be produced even 
by a pressureless and irrotational perfect fluid. In fact, the curl of the quantity pv a is still non-vanishing, even if v a 
is derived from a scalar potential. The solution of the inhomogeneous gravitational-wave equation Eq. (66) is given 
in Appendix A. 



A. Comparison with the quadrupole radiation 

We want to show how the gravitational wave source TZ^p in Eq. (68) includes the contribution by the reduced 
quadrupole moment [21] of the matter distribution expressed via comoving coordinates 



Q a p = J <Pxp(x a Xf,-^x v x v 6 a p > j . 



(69) 



First of all, let us choose the origin of our coordinates inside the mass-energy distribution described by the stress- 
energy tensor 8T l j . Let x be the vector from to the observation point P and x the vector from to the volume 
element d 3 x. On scales well inside the Hubble horizon, the solution of Eq. (66), augmented by an outgoing- wave 
boundary condition, is 

h ,fo,x) = — 4 V v % J d x |x _^ , (70) 

where the transvcrse-traceless operator is g = V a v V^p — \V a pV^ u , with V a p given by Eq. (13) and TZ" S3 by 
Eq. (68). The subscript "ret" means the quantity is to be evaluated at the retarded space-time point (rj — |x — x|/c, x). 

Our purpose is to evaluate h a p in the wave-zone, that is far outside the source region: |x| = r » x|, thus we 
expand the retarded integral Eq. (70) in powers of x/r and take only the first term of the multipole expansion 



AG 1 



a 3 J d s xn» s ^ 



(71) 



where for radially travelling waves V a p — S a p — x a xp/r 2 . Eq. (71) expresses the gravitational waves h a p in terms of 
integrals over the "stress distribution" TZ^p, while Eq. (69) represents an integral over the source "energy distribu- 
tion" . In order to make the comparison between these two equations, we need to convert the spatial components T a p 
of the stress-energy tensor in terms of the time components by means of the conservation equations = 0. Since 
Q a p in Eq. (69) is the Newtonian quadrupole and the dynamics of the tensor source is also Newtonian, we only need 
the continuity and Euler Eq. (63)- (64), by which, after some mathematical manipulations, we obtain 



/ 



d 3 x(pv a Vp) = \-^J d3S: P (x a V + XpV a ) + 2H J d 3 xp{x a Vp + XpV a ) + i / d 3 .rp (.?"<) ,s - r,i)"-f) ■ (~' 2 ') 
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and 



J d 3 xp(x a vp + Xf)V a ) = y d 3 x{px a xp) +m J d 3 x(px a x f3 ) 



(73) 



where we have dropped surface terms at infinity. 

By substituting Eq. (73) into Eq. (72), we finally find 



J d 3 x(pv a v p ) = J d 3 x (px a x p ) + ^^J d 3 xH{px a x f3 ) + 2'H J d 3 x p (x a v p + Z p v a ) 



-- / tPxp(x a dp<p + X(,ff><p) ■ 



(74) 



After substituting Eq. (74) into Eq. (68) and using again the continuity equation, in the wave zone the gravitational 
wave h a f3, to leading order in powers of 1/c and x/r, reads 



« 



d 2 Q l 



I 



drj 2 drj 

,_ / 2 
d x pi x u d ll Lp + x^d v ip - -x a d a ip 5" M 



+ lU^-^- + (m' + 12H 2 ) Q\ 



1 



2irGa 2 



d A x ( 



(75) 



ret 



Let us observe that the first line of Eq. (75) recovers the known expression of the quadrupolc radiation in the 
limit of a flat and static Universe [21, 26], while contributions on the second line derive from the back-reaction of the 
gravitational potential ip which can act as a source of gravitational waves. Moreover, on scales much smaller than the 
Hubble horizon, the last two terms on the first line in Eq. (75) can be neglected in comparison to the first one. In fact, 

the typical free fall-time of a mass distribution is proportional to p^ 1 ^ 2 , while the Hubble time goes as p^ 1 ^ 2 ', where 
"6" stands for background; thus, on small scales, where the density contrasts can be very high, the characteristic rate 
of the structure collapse is much larger than the expansion rate. This allows to drop the contributions proportional 
to H in Eq. (75) and recover the results expected well inside the horizon. 



VI. CONCLUDING REMARKS 



The main result of this paper is represented by the set of equations (25), (26), (28)-(30) and (68) , expressing metric 
perturbation in terms of the gravitational field <p, where the matter density p and the peculiar velocity w Q , satisfy 
Eqs. (24), (32) and (33). These equations, when applied in a cosmological setting characterized by a pressureless and 
irrotational fluid and a cosmological constant, provide a unified description of cosmological perturbations during their 
evolution from the linear to the highly non-linear regime. On large scales, these equations reduce to the equations 
of the first and second-order perturbation theory developed in the Poisson gauge, while, on very small scales, where 
the perturbative approach is no longer applicable, they describe the evolution of cosmological perturbations by a PN 
approximation. Indeed, we calculate the (0-0) and (0-a) components of the metric (1) up to the 1PN order, the (a-/3) 
scalar-type component up to the 2PN order, while we find for the (a-0) tensor-type component the leading-order 
source terms in powers of 1 / c. 

We also derive the generalization of the standard Euler-Poisson system of equations of Newtonian hydrodynamics, 
that consistently accounts for all the effects up to order 1/c 2 . The curl term and anisotropic stress, that produce 
vectors and tensor metric perturbations, arise already at the second perturbative order and at the strongly non-linear 
level they are dominated by the contribution of the high-density contrast and the high peculiar velocity typical of 
small-scale structures. It can be worth to stress that the quantities which source vector and tensor modes are of 
Newtonian origin on all scales, in the sense that they involve only terms that satisfy the Newtonian Euler-Poisson 
system. This result is of extreme importance in view of a possible numerical implementation of our set of equations, 
as it implies that one can compute directly vector and tensor modes starting from the outputs of N-body simulations. 

Finally, it should be stressed that our new set of equations has many possible cosmological applications such as, 
for example, the evaluation of the stochastic gravitational-wave backgrounds produced by CDM halos [27-29] and 
substructures within halos. It can be also used to improve the estimate of gravitational lensing effects and gravity- 
induced secondary CMB temperature/polarization anisotropies generated by small-scale structures [7-9]. 



12 



Acknowledgments 

We warmly thank Carlo Baccigalupi and Marco Bruni for helpful discussions. Moreover we wish to thank the 
anonymous referees which contributed with their suggestions to improve the analytic form and physical interpretation 
of our equations. 

APPENDIX A: SOLUTION OF THE INHOMOGENEOUS GRAVITATIONAL- WAVE EQUATION 

We can formally write Eq. (66) as 



(Al) 



and solve it in the two cases of matter-domination and A-domination, when the equation of state of the background 
fluid is p = wp with w = 0,-1, respectively. In order to obtain these solutions, we have Fourier-expanded the 
functions h a p(r],x) and S a p(r], x) and decomposed them in the so-called a — + , x polarization modes, as follows 



d 3 k 



h a p( V ,k)e 



ik-x 



(2tt) 3 

h a p { V , k) = h + ( V , k) p+ a p (k) + h x ( V , k) p x a p (k) , 
d 3 k 



j(?7,x) = J 



^,(,,k)e- 



S a ( V , k) = S + ( V , k) p + %(k) + S x (77, k) p* a f3(k) , 



where p +a fj(k) and p x a p{k) are the polarization tensors. 
After imposing the initial conditions 



K{Veq,ty = KiVeq,*-) = if W = , 

/i CT (r? A ,k) = /i CT ,A(k), /i' CT (?7A,k) = h' A (k) if w = -1 



(A2) 
(A3) 
(A4) 
(A5) 



(A6) 
(A7) 



(?7eg is the conformal time at matter-radiation equality, while rj A corresponds to the time when p = p A ), we obtain 
for w = 



M?? ,k) = ^% 3 



kv J Vcq kf) kq J Vet kfj 

where rj represents the present value of the conformal time. For w — — 1 we have 

M»7,k) = (v k f.h(v k ) h A,, 

- (r)k) 2 ni(r]k) h A , 



(A8) 



2ni(r7 A fc) + 77 A ni(?7Afc) k 2 /"' ni(fe»?) . 



VA 



ZjiiVAk) + VAj'AvAk) k 2 f v ji{kfj) 



VA 



K 2 


f 


H 2 k 




K 2 


r 


H$k u 





(A9) 



where /i CT ,A(k) is obtained from Eq. (A8) at r\ = ij A , while ji(r]k) and ni(r]k) are respectively the spherical Bessel and 
Neumann functions of order n = 1. 
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